Import data and define variables:

physeq = readRDS("~/Projects/ETH/Alessia/mobilome/human/human_resistome_phyloseq.rds")
out_pptx = "~/Projects/ETH/Alessia/mobilome/human/graph.pptx"

physeq %>% 
  sample_data() %>% 
  data.frame() %>% 
  mutate_at(vars(Day_of_Treatment), ~ replace(., is.na(.), -1)) -> sample_data(physeq)


physeq %>% 
  tax_table() %>% 
  data.frame() -> tax_mapping

physeq_tmp <- physeq
tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Drug_Class_multi", "Best_Hit_ARO")]
# tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Best_Hit_ARO", "Drug_Class")]

physeq_tmp %>% 
  tax_glom(., taxrank = "Best_Hit_ARO") %>% 
  tax_glom(., taxrank = "Drug_Class_multi") -> physeq_drug_multi

taxa_names(physeq_drug_multi) <- tax_table(physeq_drug_multi)[,"Drug_Class_multi"]

physeq_drug_multi

## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 19 taxa and 49 samples ]:
## sample_data() Sample Data:        [ 49 samples by 60 sample variables ]:
## tax_table()   Taxonomy Table:     [ 19 taxa by 2 taxonomic ranks ]:
## taxa are rows

physeq_tmp <- physeq
tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Drug_Class", "Best_Hit_ARO")]
# tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Best_Hit_ARO", "Drug_Class")]

physeq_tmp %>% 
  tax_glom(., taxrank = "Best_Hit_ARO") %>% 
  tax_glom(., taxrank = "Drug_Class") -> physeq_drug

taxa_names(physeq_drug) <- tax_table(physeq_drug)[,"Drug_Class"]

physeq_drug

## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 38 taxa and 49 samples ]:
## sample_data() Sample Data:        [ 49 samples by 60 sample variables ]:
## tax_table()   Taxonomy Table:     [ 38 taxa by 2 taxonomic ranks ]:
## taxa are rows

physeq_tmp <- physeq

tax_table(physeq_tmp) <- tax_table(physeq_tmp)[,c("Best_Hit_ARO", "Drug_Class")]

physeq_tmp %>% 
  tax_glom(., taxrank = "Best_Hit_ARO") -> physeq_AMRgn

taxa_names(physeq_AMRgn) <- tax_table(physeq_AMRgn)[,"Best_Hit_ARO"]

physeq_AMRgn

## phyloseq-class experiment-level object
## otu_table()   OTU Table:          [ 104 taxa and 49 samples ]:
## sample_data() Sample Data:        [ 49 samples by 60 sample variables ]:
## tax_table()   Taxonomy Table:     [ 104 taxa by 2 taxonomic ranks ]:
## taxa are rows

#define colors for Drug classes:
set.seed(123)

col <- distinctColorPalette(ntaxa(physeq_drug_multi))

names(col) <- taxa_names(physeq_drug_multi)

physeq_AMRgn %>%
  microbiome::transform(transform = "log10p") %>% 
  speedyseq::psmelt() %>% 
  left_join(tax_mapping,
            by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
            suffix = c("_x", "")) %>% 
  mutate(Abundance = na_if(Abundance, 0)) %>%
  mutate(logTPM = log10(Abundance + 1)) %>%
  separate_rows(., Drug_Class,sep = '; ',convert = TRUE) %>% 
  ggplot(mapping = aes(x = as.factor(Day_of_Treatment),
                       y = Best_Hit_ARO ,
                       fill = logTPM,
  )) +
  facet_grid(Drug_Class ~ Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free") + #,  labeller = labeller(Drug_Class_multi = label_wrap_gen(width = 15)))  +
  # facet_grid(. ~ Reactor_Treatment, scales = "free", space = "free", drop = TRUE,labeller = labeller(Drug.Class = label_wrap_gen(width = 15))) +
  #we adjust the vertical labels such that it is printed correctly. 
  geom_tile() +
  theme_light()  +
  # scale_fill_viridis_c(name = "TPM",
  #                      na.value = "transparent", trans = scales::pseudo_log_trans(sigma = 0.1),
  #                      # trans =  scales::pseudo_log_trans(sigma = 1),
  #                      breaks = c(1, 10, 50, 100, 400), labels = c(1, 10, 50, 100, 400),
  #                      limits = c(0, 500)) + #need to make sure that we include max value of tpm as upper limit
  scale_fill_viridis_c(na.value = "transparent") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        strip.text.x = element_text(size=8),
        #change the size of the label and the angle so the text is printed horizontally
        strip.text.y = element_text(size = 7,angle=0)) -> perma_bbuble

perma_bbuble +
  theme(legend.position = "bottom") -> perma_bbuble 

perma_bbuble +
  scale_y_discrete(label = function(x) stringr::str_trunc(x, 34, side = "center")) ->perma_bbuble

perma_bbuble

perma_bbuble %>%
  export::graph2ppt(append = TRUE, width = 20, height = 50,
                    file = out_pptx)

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

physeq_AMRgn %>%
  microbiome::transform(transform = "log10p") %>% 
  speedyseq::psmelt() %>% 
  left_join(tax_mapping,
            by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
            suffix = c("_x", "")) %>% 
  mutate(Best_Hit_ARO = fct_reorder(Best_Hit_ARO, Abundance, .desc=FALSE)) %>% 
  mutate(Abundance = na_if(Abundance, 0)) %>%
  mutate(logTPM = log10(Abundance + 1))  %>% 
  # df$Item = factor(df$Item, levels = a_lot$Item[order(a_lot$percent)])

# tmp$variable <- factor(melted$variable, levels=rev(levels(melted$variable)))
  # mutate(Best_Hit_ARO =  fct_reorder(Best_Hit_ARO, logTPM)) %>% 
  # mutate(Best_Hit_ARO =  fct_relevel(Best_Hit_ARO, logTPM)) %>% 
  # separate_rows(., Drug_Class,sep = '; ',convert = TRUE) %>% 
  ggplot(mapping = aes(x = as.factor(Day_of_Treatment),
                       y = Best_Hit_ARO,
                       fill = logTPM,
  )) +
  facet_grid(.  ~  Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free") + #,  labeller = labeller(Drug_Class_multi = label_wrap_gen(width = 15)))  +
  # facet_grid(. ~ Reactor_Treatment, scales = "free", space = "free", drop = TRUE,labeller = labeller(Drug.Class = label_wrap_gen(width = 15))) +
  #we adjust the vertical labels such that it is printed correctly. 
  geom_tile() +
  theme_light()  +
  # scale_fill_viridis_c(name = "TPM",
  #                      na.value = "transparent", trans = scales::pseudo_log_trans(sigma = 0.1),
  #                      # trans =  scales::pseudo_log_trans(sigma = 1),
  #                      breaks = c(1, 10, 50, 100, 400), labels = c(1, 10, 50, 100, 400),
  #                      limits = c(0, 500)) + #need to make sure that we include max value of tpm as upper limit
  scale_fill_viridis_c(na.value = "transparent") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        strip.text.x = element_text(size=8),
        #change the size of the label and the angle so the text is printed horizontally
        strip.text.y = element_text(size = 7,angle=0)) -> perma_bbuble

perma_bbuble +
  theme(legend.position = "bottom") +
  scale_y_discrete(label = function(x) stringr::str_trunc(x, 34, side = "center"))  -> perma_bbuble

perma_bbuble

perma_bbuble %>%
  export::graph2ppt(append = TRUE, width = 20, height = 28,
                    file = out_pptx)

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

physeq_AMRgn %>%
  microbiome::transform(transform = "log10p") %>% 
  speedyseq::psmelt() %>% 
  left_join(tax_mapping,
            by = c("Best_Hit_ARO" = "Best_Hit_ARO"),
            suffix = c("_x", "")) %>% 
  mutate(Best_Hit_ARO = fct_reorder(Best_Hit_ARO, Abundance, .desc=FALSE)) %>%
  mutate(Abundance = na_if(Abundance, 0)) %>%
  mutate(logTPM = log10(Abundance + 1))  %>% 
  # df$Item = factor(df$Item, levels = a_lot$Item[order(a_lot$percent)])

# tmp$variable <- factor(melted$variable, levels=rev(levels(melted$variable)))
  # mutate(Best_Hit_ARO =  fct_reorder(Best_Hit_ARO, logTPM)) %>% 
  # mutate(Best_Hit_ARO =  fct_relevel(Best_Hit_ARO, logTPM)) %>% 
  # separate_rows(., Drug_Class,sep = '; ',convert = TRUE) %>% 
  ggplot(mapping = aes(x = as.factor(Day_of_Treatment),
                       y = Best_Hit_ARO,
                       fill = logTPM,
  )) +
  facet_grid(Drug_Class_multi  ~  Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free") + #,  labeller = labeller(Drug_Class_multi = label_wrap_gen(width = 15)))  +
  # facet_grid(. ~ Reactor_Treatment, scales = "free", space = "free", drop = TRUE,labeller = labeller(Drug.Class = label_wrap_gen(width = 15))) +
  #we adjust the vertical labels such that it is printed correctly. 
  geom_tile() +
  theme_light()  +
  # scale_fill_viridis_c(name = "TPM",
  #                      na.value = "transparent", trans = scales::pseudo_log_trans(sigma = 0.1),
  #                      # trans =  scales::pseudo_log_trans(sigma = 1),
  #                      breaks = c(1, 10, 50, 100, 400), labels = c(1, 10, 50, 100, 400),
  #                      limits = c(0, 500)) + #need to make sure that we include max value of tpm as upper limit
  scale_fill_viridis_c(na.value = "transparent") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        strip.text.x = element_text(size=8),
        #change the size of the label and the angle so the text is printed horizontally
        strip.text.y = element_text(size = 7,angle=0)) -> perma_bbuble

perma_bbuble +
  theme(legend.position = "bottom") +
  scale_y_discrete(label = function(x) stringr::str_trunc(x, 34, side = "center"))  -> perma_bbuble

perma_bbuble

perma_bbuble %>%
  export::graph2ppt(append = TRUE, width = 20, height = 28,
                    file = out_pptx)

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

# see Sneha's & Hannah's code for "plot.df.sep <- separate_rows(plot.df, Drug.Class,sep = '; ',convert = TRUE)"

physeq %>%
  microbiome::transform(transform = "log10p") %>% 
  speedyseq::psmelt() %>% 
  separate_rows(., Resistance_Mechanism,sep = '; ',convert = TRUE) %>% 
  ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Resistance_Mechanism)) +   
  geom_bar( stat = "identity", colour="black", size=0.05) +
  facet_grid(Drug_Class_multi  ~  Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
  ggtitle("ARG Type Abundance per Sample") +
  xlab("Day (Treatment)")  +
  ylab("log10 TPM") + 
  theme_light() -> p_mec

p_mec %>% 
  ggpubr::set_palette("jco") -> p_mec

p_mec

p_mec %>%
  export::graph2ppt(append = TRUE, width = 20, height = 18,
                    file = out_pptx)

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

physeq_drug_multi %>%
  microbiome::transform(transform = "log10p") %>% 
  speedyseq::psmelt() %>% 
  ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Drug_Class_multi)) +   
  geom_bar( stat = "identity", colour="black", size=0.05) +
  scale_fill_manual(values=col) +
  #geom_bar(stat="identity",colour=NA,size=0) +
  facet_grid( ~  Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
  ggtitle("ARG Type Abundance per Sample") +
  xlab("Day (Treatment)")  +
  ylab("log10 TPM") + 
  theme_light() -> p_drug

p_drug

p_drug %>%
  export::graph2ppt(append = TRUE, width = 18, height = 9,
                    file = out_pptx)

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

physeq_drug_multi %>%
  microbiome::transform(transform = "log10p") %>% 
  speedyseq::psmelt() %>% 
  ggplot(aes(x = as.factor(Day_of_Treatment), y = Abundance, fill = Drug_Class_multi)) +   
  geom_bar( stat = "identity", position="fill", colour="black", size=0.05) +
  scale_fill_manual(values=col) +
  #geom_bar(stat="identity",colour=NA,size=0) +
  facet_grid( ~  Treatment + Fermentation + Antibiotic_mg.mL , scales = "free", space = "free_x") +
  ggtitle("ARG Type Abundance per Sample") +
  xlab("Day (Treatment)")  +
  ylab("Proportion") +
  theme_light() -> p_drug_prop

p_drug_prop

p_drug_prop %>%
  export::graph2ppt(append = TRUE, width = 18, height = 9,
                    file = out_pptx)

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

ARG Alpha Diversity

# extract ARG abundance info from phyloseq object
physeq_AMRgn %>%
  otu_table() %>%
  as.data.frame() %>%
  t() -> resistome_copy

# extract sample info
physeq_AMRgn %>%
  sample_data() %>%
  as.data.frame() -> treatment_info

# add sample info to resistome abundance data and get richeness
resistome_copy %>%
  specnumber() %>%
  cbind(., treatment_info) %>% 
  mutate(metric = "Richness") %>% 
  rownames_to_column('sample_id') -> ARG_richeness #%>%
# filter(., Treatment != "DONOR")-> ARG_richeness

colnames(ARG_richeness)[2] <- "value"

# resistome_copy %>%
#   diversity() %>%
#   cbind(., treatment_info) %>% 
#   mutate(metric = "Shannon") %>% 
#   rownames_to_column('sample_id') -> ARG_shannon #%>%
# 
# colnames(ARG_shannon)[2] <- "value"
# 
# ARG_shannon %>% 
#   rbind(ARG_richeness) -> alpha_resistome

ARG_richeness -> alpha_resistome

# levels(alpha_resistome$Day_of_Treatment)[alpha_resistome$Day_of_Treatment] %>%
#   as.numeric() -> alpha_resistome$Day_of_Treatment



# plot
ggplot(alpha_resistome, aes(x = Reactor_Treatment, y = value)) + 
  geom_boxplot(aes(color=Reactor_Treatment),
               outlier.shape = NA,
               outlier.colour = NA,
               outlier.size = 0) + 
  geom_jitter(aes(color=Reactor_Treatment), alpha = 0.4) +
  theme(legend.position="none") +
  xlab("Treatment") +
  ylab("ARG Richness") +
  theme_light() + scale_color_viridis_d(na.value = "black") +
  facet_grid(metric ~ .)

# plot panels for each treatment
alpha_resistome %>% 
  # filter(!is.na(Fermentation)) %>% 
  # filter(Treatment != "DONOR") %>%
  ggplot(aes(x = Day_of_Treatment, y = value, color = Treatment)) +
  geom_point() + 
  geom_line(linetype = "dashed",  size = 0.25) +
  labs(title = "AMR Gene Richness",
       y = "Richness", x = "Days of Treatment") +
  scale_color_viridis_d(na.value = "black") +
  facet_grid(.  ~  Treatment + Fermentation + Antibiotic_mg.mL , scales = "fixed", space = "free_x") +
  theme_light() -> p_alpha

p_alpha

## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

p_alpha %>% 
  export::graph2ppt(append = TRUE, width = 18, height = 6,
                    file = out_pptx)

## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

Beta Diversity:

# sample_data(physeq)$Reactor_Treatment <- fct_relevel(sample_data(physeq)$Reactor_Treatment, "DONOR", "CR_UNTREATED", "TR1_CTX+HV292.1",
#                                                      "TR2_CTX","TR3_HV292.1", "TR4_VAN", "TR5_VAN+CCUG59168", "TR6_CCUG59168") 
# 
# sample_data(physeq)$Treatment <- fct_relevel(sample_data(physeq)$Treatment, "DONOR", "UNTREATED",  "CTX+HV292.1", "CTX","HV292.1","VAN+CCUG59168", "VAN",  "CCUG59168") 


# physeq %>% 
#   rarefy_even_depth(sample.size = 43,
#                     rngseed = 123) -> phyloseq_rare

physeq_AMRgn %>%
  phyloseq_compute_bdiv(norm = "pc",
                        phylo = FALSE,
                        seed = 123) -> bdiv_list

## Loading required package: ape

## Loading required package: GUniFrac

## Registered S3 method overwritten by 'rmutil':
##   method         from
##   print.response httr

physeq_AMRgn  %>%
  # subset_samples(Treatment != "DONOR") %>%
  phyloseq_plot_bdiv(dlist = bdiv_list,
                     seed = 123,
                     axis1 = 1,
                     axis2 = 2)  -> pcoa

## [1] "bray"
## [1] "sorensen"
## [1] "bjaccard"
## [1] "wjaccard"

# phyloseq_plot_bdiv(bdiv_list,
#                    # m = "CoDa",
#                    seed = 123) -> coda
# 
pcoa$wjaccard$layers = NULL

pcoa$wjaccard + geom_point(size=3,
                           aes(color = Treatment, 
                               fill = NULL,
                               shape = Fermentation %>%  as.factor(),
                               alpha = Day_of_Treatment)) + 
  theme_light() +
  geom_path(arrow = arrow(type = "open", angle = 30, length = unit(0.15, "inches")),
            size = 0.05, linetype = "dashed", inherit.aes = TRUE, aes(group=Reactor_Treatment), show.legend = FALSE) +
  scale_alpha_continuous(range=c( 0.9, 0.3)) +
  scale_color_viridis_d(na.value = "black") + 
  scale_fill_viridis_d(na.value = "black") + 
  scale_shape_manual(values = c(15, 19), na.value =  17) + 
  theme_classic() -> p1 # +
  # labs(col=NULL, fill = NULL, shape = NULL) + guides(shape=TRUE) -> p1

p1

p1 %>% 
  export::graph2ppt(append = TRUE, width = 8, height = 6,
                    file = out_pptx)

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

physeq_AMRgn %>% 
  tax_table() %>% 
  data.frame() %>% 
  rownames_to_column("id") %>% 
  mutate(gene = id) %>% 
  column_to_rownames("id") %>% 
  as.matrix %>% 
  tax_table() -> tax_table(physeq_AMRgn)



physeq_AMRgn  %>%
  # subset_samples(Treatment != "DONOR") %>% 
  phyloseq_add_taxa_vector(dist = bdiv_list$wjaccard,
                           phyloseq = .,
                           figure_ord = p1,
                           tax_rank_plot = "Best_Hit_ARO", taxrank_glom = "Best_Hit_ARO",
                           top_r = 10, fact = 0.6) -> pco_env


pco_env$plot

## Warning: Removed 3 rows containing missing values (geom_text_repel).

pco_env$plot %>% 
  export::graph2ppt(append = TRUE, width = 8, height = 6,
                    file = out_pptx)

## Warning: Removed 3 rows containing missing values (geom_text_repel).

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

pco_env$signenvfit %>% 
  DT::datatable()

physeq_drug_multi  %>%
  # subset_samples(Treatment != "DONOR") %>% 
  phyloseq_add_taxa_vector(dist = bdiv_list$wjaccard,
                           phyloseq = .,
                           figure_ord = p1,
                           tax_rank_plot = "Drug_Class_multi", taxrank_glom = "Drug_Class_multi",
                           top_r = 10, fact = 0.6) -> pco_env


pco_env$plot

pco_env$plot %>% 
  export::graph2ppt(append = TRUE, width = 8, height = 6,
                    file = out_pptx)

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

pco_env$signenvfit %>% 
  DT::datatable()

Make the heatmap pretty

physeq %>% 
  # speedyseq::tax_glom("Best_Hit_ARO") %>% 
  speedyseq::psmelt() %>% 
  rename(TPM = Abundance) %>% 
  mutate(TPM = na_if(TPM, 0)) %>%
  # mutate(logTPM = log10(TPM + 1)) %>% 
  # filter(Drug_Class %in% c("cephalosporin", "glycopeptide")) %>% 
  # mutate(Drug_Class = fct_relevel(Drug_Class, c("cephalosporin", "glycopeptide"))) %>% 
  filter(Best_Hit_ARO %in% c("vanZA", "vanYA", "vanXA", "vanA", "vanHA", "vanRA", "CTX-M-1")) %>% 
  # mutate(Classification = fct_relevel(Classification, c("MAG", "Chromosome", "Plasmid", "Uncertain - plasmid or chromosomal"))) %>% 
  ggplot(mapping = aes(x = Day_of_Treatment ,
                       y = Best_Hit_ARO ,
                       fill = TPM,
  )) +
  facet_grid(Drug_Class ~  Treatment + Fermentation + Antibiotic_mg.mL, scales = "free_y",  space = "free_y", drop = TRUE) +
  # facet_grid(. ~ Treatment,  scales = "free", space = "free", drop = TRUE) +
  # facet_grid(. ~ Reactor_Treatment, scales = "free", space = "free", drop = TRUE,labeller = labeller(Drug_Class = label_wrap_gen(width = 15))) +
  #we adjust the vertical labels such that it is printed correctly. 
  geom_tile() +
  theme_bw()  +
  # scale_fill_viridis_c(name = "TPM",
  #                      na.value = "transparent", trans = scales::pseudo_log_trans(sigma = 0.1),
  #                      # trans =  scales::pseudo_log_trans(sigma = 1),
  #                      breaks = c(1, 10, 50, 100, 400), labels = c(1, 10, 50, 100, 400),
  #                      limits = c(0, 500)) + #need to make sure that we include max value of tpm as upper limit
  scale_fill_viridis_c(na.value = "transparent") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        strip.text.x = element_text(size=8),
        #change the size of the label and the angle so the text is printed horizontally
        strip.text.y = element_text(size = 7,angle=0)) -> perma_bbuble

perma_bbuble +
  theme(legend.position = "bottom") -> perma_bbuble

perma_bbuble

perma_bbuble %>% 
  export::graph2ppt(append = TRUE, width = 10, height = 4,
                    file = out_pptx)

## Exported graph as ~/Projects/ETH/Alessia/mobilome/human/graph.pptx

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] GUniFrac_1.4         ape_5.6              gdtools_0.2.3       
##  [4] reshape2_1.4.4       scales_1.1.1         randomcoloR_1.1.0.1 
##  [7] vegan_2.5-7          lattice_0.20-45      permute_0.9-5       
## [10] RColorBrewer_1.1-2   microViz_0.9.0       here_1.0.1          
## [13] ggrepel_0.9.1        speedyseq_0.5.3.9018 phyloseq_1.36.0     
## [16] forcats_0.5.1        stringr_1.4.0        dplyr_1.0.7         
## [19] purrr_0.3.4          readr_2.1.0          tidyr_1.1.4         
## [22] tibble_3.1.6         ggplot2_3.3.5        tidyverse_1.3.1.9000
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1           uuid_1.0-3             backports_1.4.1       
##   [4] systemfonts_1.0.2      plyr_1.8.6             igraph_1.2.10         
##   [7] splines_4.1.2          crosstalk_1.2.0        GenomeInfoDb_1.28.4   
##  [10] digest_0.6.29          foreach_1.5.1          htmltools_0.5.2       
##  [13] lmerTest_3.1-3         fansi_0.5.0            magrittr_2.0.1        
##  [16] cluster_2.1.2          tzdb_0.2.0             openxlsx_4.2.4        
##  [19] Biostrings_2.60.2      modelr_0.1.8           matrixStats_0.61.0    
##  [22] stabledist_0.7-1       officer_0.4.0          colorspace_2.0-2      
##  [25] rvest_1.0.2            haven_2.4.3            xfun_0.28             
##  [28] crayon_1.4.2           RCurl_1.98-1.5         jsonlite_1.7.2        
##  [31] lme4_1.1-27.1          survival_3.2-13        iterators_1.0.13      
##  [34] glue_1.6.0             rvg_0.2.5              gtable_0.3.0          
##  [37] zlibbioc_1.38.0        XVector_0.32.0         V8_3.6.0              
##  [40] car_3.0-11             Rhdf5lib_1.14.2        BiocGenerics_0.38.0   
##  [43] abind_1.4-5            DBI_1.1.1              rstatix_0.7.0         
##  [46] Rcpp_1.0.7             viridisLite_0.4.0      xtable_1.8-4          
##  [49] clue_0.3-60            foreign_0.8-81         DT_0.20               
##  [52] stats4_4.1.2           timeSeries_3062.100    htmlwidgets_1.5.4     
##  [55] httr_1.4.2             ellipsis_0.3.2         spatial_7.3-14        
##  [58] pkgconfig_2.0.3        farver_2.1.0           sass_0.4.0            
##  [61] dbplyr_2.1.1           utf8_1.2.2             tidyselect_1.1.1      
##  [64] labeling_0.4.2         rlang_0.4.12           munsell_0.5.0         
##  [67] cellranger_1.1.0       tools_4.1.2            cli_3.1.0             
##  [70] generics_0.1.1         devEMF_4.0-2           ade4_1.7-18           
##  [73] export_0.3.0           broom_0.7.11           evaluate_0.14         
##  [76] biomformat_1.20.0      fastmap_1.1.0          yaml_2.2.1            
##  [79] knitr_1.36             fs_1.5.2               zip_2.2.0             
##  [82] rgl_0.107.14           nlme_3.1-153           xml2_1.3.2            
##  [85] compiler_4.1.2         rstudioapi_0.13        curl_4.3.2            
##  [88] ggsignif_0.6.3         reprex_2.0.1           statmod_1.4.36        
##  [91] statip_0.2.3           bslib_0.3.1            stringi_1.7.6         
##  [94] highr_0.9              modeest_2.4.0          stargazer_5.2.2       
##  [97] fBasics_3042.89.1      Matrix_1.3-4           nloptr_1.2.2.2        
## [100] ggsci_2.9              microbiome_1.14.0      multtest_2.48.0       
## [103] vctrs_0.3.8            pillar_1.6.4           lifecycle_1.0.1       
## [106] rhdf5filters_1.4.0     jquerylib_0.1.4        data.table_1.14.2     
## [109] bitops_1.0-7           flextable_0.6.9        stable_1.1.4          
## [112] R6_2.5.1               rio_0.5.27             IRanges_2.26.0        
## [115] codetools_0.2-18       boot_1.3-28            MASS_7.3-54           
## [118] assertthat_0.2.1       rhdf5_2.36.0           rprojroot_2.0.2       
## [121] withr_2.4.3            S4Vectors_0.30.2       GenomeInfoDbData_1.2.6
## [124] mgcv_1.8-38            parallel_4.1.2         hms_1.1.1             
## [127] rpart_4.1-15           timeDate_3043.102      grid_4.1.2            
## [130] minqa_1.2.4            rmarkdown_2.11         rmutil_1.1.5          
## [133] carData_3.0-4          Rtsne_0.15             ggpubr_0.4.0          
## [136] numDeriv_2016.8-1.1    Biobase_2.52.0         lubridate_1.8.0       
## [139] base64enc_0.1-3

Resistome_visualization

Hannah Li Hägi & FC

January 07, 2022

Load required packages

Source required functions

Import data and define variables: